# -*- coding: utf-8 -*-
"""
Created on 2023-07-20
势流理论圆柱绕流气动力分析 - 带环量修正
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import simps
import pandas as pd
class AerodynamicIntegrator:
    @staticmethod
    def validate_dAlembert_paradox(theta, Cp, tolerance=1e-6):
        """
        验证达朗贝尔佯谬（阻力是否为零）
        参数：
            theta : 角度数组
            Cp    : 压力系数数组
            tolerance : 允许的数值误差
        返回：
            is_paradox_hold : 是否满足佯谬（True/False）
            CD : 计算得到的阻力系数
        """
        integrand_drag = Cp * np.cos(theta)
        CD = (1/(2*np.pi)) * simps(integrand_drag, theta)
        is_paradox_hold = abs(CD) < tolerance
        return is_paradox_hold, CD
        if __name__ == "__main__":
    params = FlowParameters()
    calculator = PotentialFlowCalculator()
    integrator = AerodynamicIntegrator()
    
    # 验证无环量时的达朗贝尔佯谬
    Cp0 = calculator.pressure_coefficient(params.theta, Gamma=0, params=params)
    is_paradox, CD = integrator.validate_dAlembert_paradox(params.theta, Cp0)
    
    print("\n达朗贝尔佯谬验证结果：")
    print(f"理论预测阻力 = 0 | 数值计算阻力 = {CD:.2e}")
    print(f"是否满足佯谬？ {'是' if is_paradox else '否'}")
    
    # 可视化压力分布对称性
    plt.figure(figsize=(8, 4))
    plt.plot(params.theta, Cp0 * np.cos(params.theta), 'b-', label=r'$C_p \cos\theta$ (阻力积分项)')
    plt.axhline(0, color='k', linestyle='--')
    plt.title("达朗贝尔佯谬：阻力积分项的对称性")
    plt.xlabel(r"$\theta$ (rad)")
    plt.legend()
    plt.grid()
    plt.show()
# ======================== 物理模型参数类 ========================
class FlowParameters:
    """定义流体力学计算参数"""
    def __init__(self):
        self.a = 1.0          # 圆柱半径 (m)
        self.U = 1.0          # 来流速度 (m/s)
        self.rho = 1.225      # 空气密度 (kg/m³)
        self.Gamma_list = [    # 环量参数 (m²/s)
            0, 
            2*np.pi*self.a*self.U, 
            4*np.pi*self.a*self.U
        ]
        self.theta = np.linspace(0, 2*np.pi, 1000)  # 角度离散化

# ======================== 理论计算模块 ========================
class PotentialFlowCalculator:
    """势流理论气动力计算器"""
    @staticmethod
    def pressure_coefficient(theta, Gamma=None, params=None):
        """
        计算压力系数分布
        参数：
            theta : 方位角数组 (rad)
            Gamma : 环量值 (m²/s)
            params: FlowParameters实例
        返回：
            Cp : 压力系数数组
        """
        a, U = params.a, params.U
        
        if Gamma is None or Gamma == 0:
            # 无环量情况
            return 1 - 4*np.sin(theta)**2
        else:
            # 含环量修正
            term1 = -4*np.sin(theta)**2
            term2 = -(2*Gamma/(np.pi*a*U)) * np.sin(theta)
            term3 = -(Gamma/(2*np.pi*a*U))**2
            return 1 + term1 + term2 + term3

    @staticmethod
    def theoretical_lift(Gamma, params):
        """库塔-茹科夫斯基升力理论值"""
        return -Gamma / (params.a * params.U)

# ======================== 数值积分模块 ========================
class AerodynamicIntegrator:
    """气动力数值积分器"""
    @staticmethod
    def compute_forces(theta, Cp):
        """
        计算无量纲气动力系数
        参数：
            theta : 方位角数组 (rad)
            Cp    : 压力系数数组
        返回：
            CL : 升力系数
            CD : 阻力系数
        """
        integrand_lift = Cp * np.sin(theta)
        integrand_drag = Cp * np.cos(theta)
        
        CL = (1/(2*np.pi)) * simps(integrand_lift, theta)
        CD = (1/(2*np.pi)) * simps(integrand_drag, theta)
        return np.round(CL, 6), np.round(CD, 6)  # 保留6位小数

# ======================== 可视化模块 ========================
class ResultVisualizer:
    """结果可视化工具"""
    @staticmethod
    def plot_pressure_distribution(theta, Cp_data, Gamma_list):
        """绘制压力系数分布曲线"""
        plt.figure(figsize=(12, 6))
        for i, (Gamma, Cp) in enumerate(Cp_data):
            plt.subplot(1, len(Cp_data), i+1)
            plt.plot(theta, Cp, 'r-', linewidth=2)
            plt.title(f"Γ = {Gamma:.1f} m²/s", fontsize=12)
            plt.xlabel(r'$\theta$ (rad)', fontsize=10)
            plt.ylabel('$C_p$', fontsize=10)
            plt.grid(True, linestyle='--', alpha=0.7)
        plt.tight_layout()
        plt.show()

# ======================== 主程序 ========================
if __name__ == "__main__":
    # 初始化参数
    params = FlowParameters()
    calculator = PotentialFlowCalculator()
    integrator = AerodynamicIntegrator()
    
    # 计算所有工况
    results = []
    Cp_data = []
    for Gamma in params.Gamma_list:
        # 计算压力系数
        Cp = calculator.pressure_coefficient(params.theta, Gamma, params)
        
        # 计算气动力系数
        CL, CD = integrator.compute_forces(params.theta, Cp)
        theory_CL = calculator.theoretical_lift(Gamma, params)
        
        # 保存结果
        results.append({
            'Gamma': Gamma,
            'Theory_CL': theory_CL,
            'Numerical_CL': CL,
            'Numerical_CD': CD,
            'Relative_Error (%)': abs((CL - theory_CL)/theory_CL)*100 if theory_CL !=0 else 0
        })
        Cp_data.append( (Gamma, Cp) )
    
    # 结果表格输出
    df = pd.DataFrame(results)
    print("\n" + "="*60)
    print("气动力系数对比分析表")
    print("="*60)
    print(df.round(4).to_markdown(index=False))
    
    # 可视化
    ResultVisualizer.plot_pressure_distribution(params.theta, Cp_data, params.Gamma_list)

    # 理论局限性分析
    print("\n" + "="*60)
    print("理论模型局限性总结")
    print("="*60)
    print("1. 无粘假设导致零阻力 (数值阻力: {:.1e}量级)".format(df['Numerical_CD'].abs().max()))
    print("2. 环量需人为指定，无法反映启动涡生成机制")
    print("3. 无法预测流动分离现象 (实际Re>40会出现卡门涡街)")